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Abstract 

Combining several independent measurements of the same physical quantity is one 
of the most important tasks in metrology. Small samples, biased input estimates, 
not always adequate reported uncertainties, and unknown error distribution make a 
rigorous solution very difficult, if not impossible. For this reason, many methods to 
compute a common mean and its uncertainty were proposed, each with own advantages 
and shortcomings. Most of them are variants of the weighted average (WA) approach 
with different strategies to compute WA and its standard deviation. Median estimate 
became also increasingly popular during recent years. In this paper, these two methods 
in most widely used modifications are compared using simulated and real data. To 
overcome some problems of known approaches to compute the WA uncertainty, a new 
combined estimate has been proposed. It has been shown that the proposed method 
can help to obtain more robust and realistic estimate suitable for both consistent and 
discrepant measurements. 

1 Introduction 

Computation of a common mean (CM) of several independent measurements of a physical 
quantity is a common procedure in scientific analysis. Typical task include combining results 
obtained by different analysts or results based on different measurement methods or results 
of several series of measurements, etc. One of the most important applications required 
computation of a CM is derivation of the best estimate of physical constants. Both accurate 
CM value and its realistic uncertainty are equally important in such computations. 

Input information consists of measured values Xi, i = 1, . . . ,n with associated uncertain- 
ties Si, usually standard deviations (STD), and correlations between Xi are not available. 
From statistical point of view, it is a classical case of direct measurements of unequal preci- 
sion. Due to lack of needed information we have to treat them as uncorrelated. Combination 
procedure aimed to get an estimate x of the CM x with associated uncertainty a, which ad- 
equately reflects the scatter and uncertainties of input measurements. 

There is no unambiguous solution of this problem. It is well known the the classical 
weighted average (WA) is unbiased CM estimate with minimum variance, provided Xj are 
independent unbiased estimators, and are true variances. A solution of the problem 
is not straightforward if true variances are unknown. Numerous papers are devoted to 
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computation of a CM and its uncertainty; see e.g. iGraybill and Deall ( 1l959l ): ISinhal ( 1l985[ ): 



Witkovsky and Wimmerl (120011 ) : IZhang] (120061 ) and references therein. Unfortunately, results 
are obtained under rather strong assumptions: the Xi are unbiased estimates of x, Sj are 
normally distributed, and number of measurements rij formed each Xi is known. Evidently, 
these assumptions are hardly met in scientific data analysis. The situation becomes even 
more complicated because small samples are are the rule rather than the exception in such 
tasks, which makes it difficult to efficiently apply most of statistical methods generally used. 

For these reasons many alte r nativ e approaches , not a lways strongly justified, are mostly 
used, see e.g. iMacMahon et all (120041 ). I Chen et all (1201 ll ). and Dataplot software documen- 
tatioE0- We do not aim this paper at investigation of all the existing methods. Such a 
task seems to be impractical because some methods cannot be used in our case because we 
lack needed information, some are developed for specific app lications, and some approaches 
may be too complicated for routine use, e.g. bootstr apping (iHelend . l2007l ) or total median 
(iFigueiredo and Gomesl . l2004l : I Cox and Harris! . |200J) . 

We consider two basic methods most commonly used: WA and median. As to the former, 
several approaches to compute its uncertainty (STD) are proposed in literature. This paper 
is devoted to comparison of this methods based on simulated and real data. Besides, a new 
combined approach is proposed to c ompute the WA uncertainty. This approach was proposed 
for the first time in iMalkinI (l2001bl ) and has shown to be useful in practical applications suc h 
as computation of Earth orientation paramet er (^EOP) combined solution (iMalkinl . l2001al ). 
radio source position catalo gues cornbinati on (ISokolova and Malkinl . 120071 ) . analysis of radio 
sourc e position time series (iMalkinl . |2008|), and modeling the galactic aberration (iMalkinl . 



20111). Special attention is given to analysis of small statistical samples. 



2 Computation of CM 
2.1 Basic WA estimators 

The WA estimator is most widely used in various scientific and practical applications. Let 
we have n values Xi with ass ociated s t andard deviations Sj, i = 1 . . . n. Then we can compute 
the following statistics (e.g. iBrandtl (119991 ): iBevington and RobinsonI (120031 )) 
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In result, we have the classical WA estimate x^ with weights inversely proportional to 
variances of input measurements sf. We also can compute the following statistics also used 
as a measure of goodness of fit 
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H = "^Piixi - 



i=l 
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i=l 



1 Xi X^. 



where H has a distribution with n — 1 degree of freedom (dof) if are theoretical 
variances. In practice, rule, sample variances, but this fact is usually ignored. 



^http: / /www. itl.nist.gov/div898 / software/dataplot / 
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The H statistics is an indicator of consistency of input measurements. If the measurements 



are consistent, the value 



X 



H 



^ f 1 ^3) 

aoj n — 1 

is close to unity. Otherwise one can assume that the input measurements Xi have systematic 
errors or Si are underestimated. 

The question is how to estimate the standard error a of the mean? Two main approaches 
can be applied to compute a. The classical WA estimate is 
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Least squares approach leads to alternative estimate of the WA uncertainty. The solution 
of the least square problem Xi = x + Ei with weights Pi gives the same WA estimate x, but 
another estimate of its uncertainty: 
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This estimate of the WA uncertainty is, in fact, (Xi estimate scaled in such a way to make 



X /dof close to unity: 



H 



0-2 = (Tl 
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which is equivalent to scaling of input Sj by the factor ^H/{n — 1) as recommended by 
Rosenfeld et all ( 119671 ) and iBrandtl (Il999l ). However, such a scaling makes resulting a esti- 



mate independent of the "scale" of input variances s^, and dependent only on their ratio. 

So, both approaches give the same estimate for the WA, but different estimates for the 
WA uncertainty. The first value ai depends on Si and does not depend on the scatter of the 
input values Xj. On the other hand, a2 depends on relative values of input variances sf and 
the scatter of Xj. Difference between the two cr estimates may be attributed to systematic 
errors in Xi or underestimated Sj. 

Theoretically, choice between ai and a2 depends on whether the scatter of Xj is a result 
of random error or there are systematic differences between estimates Xj. Obviously, both 
effects are present in most of practical applications. This is a well recognized problem in data 
analysis, and its rigorous solution is hardly possible due to generally biased input estimates, 
not always adequate reported uncertainties, and unknown error distribution. 

In practice, if Xj are close each other and Si are greater than the scatter of Xj, it seems 
reasonable to use cxi. Otherwise, if Si are much less than the scatter of the input measure- 
ments, (72 estimate seems to be more adequate to the data. Indeed such a way to choose the 
best estimate for the WA uncertainty cannot be considered satisfactory. 

A possi ble prac t ical, s tatistically based approach has been proposed by Rosenfeld et al. 



(1967) and iBrandtl (119991 ). According to this approach, criteria is used to decide whether 
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the scatter of Xj is a result of random errors. First, both uncertainty estimates ai and (T2 
are computed. Then the final WA uncertainty is taken as 



ai, ifH<xHQ,n-l) 
a2, ifH>x'{Q,n-l) 



(7) 



where Q is a significance level. One can see that to a first approximation as = a\ for 
consistent measurements and (T3 = 02 for discrepant ones, and given Q value is used to 
distinguish between them. As a consequence, substantially different a estimates can be 
obtained for the same input measurer nents {xi^S j ] but specifying different Q. 

Similar approach is discussed by Bich et all ( 2002 ). It recommends accept a\ as the 
estimate of the CM mean uncertainty if consistency check H < x^((5,n — 1) passed at a 
significance level Q=5%; otherwise supplement studies should be performed, such check of 
outliers, investigation of input data, etc. Unfortunately, the latter is generally not feasible. 



2.2 Combined WA uncertainty estimator 

As pointed out in the previous section, the literature recommends using either ai or (T2 
depending on some criteria, which can lead to ambiguous results, keeping in mind that these 
two estimates may differ by several times. So, a more robust a estimate is desirable for 
practical use, which would account of both the scatter of Xi and their uncertainties Sj. After 
investigation of behavior of all three estimates using simulated and real data, and supplement 
tests we decided in favor of combined estimate computed by simple formula 

ac = ^(yl + al . (8) 

As a variant, a2 computed with unit weights can be used, which provides clear separation 
of impacts at combined estimate cxc from the uncertainties (ai) and scatter ((T2) of the input 
measurements. However, in this case, the result is generally more sensitive to measurements 
suspected to be outliers. 

Unfortunately, we cannot suggest a rigorous theoretical ground of this approach, which 
is common for other practical recommendations too. Our considerations are as follows. 
Suppose we can represent each input value as Xi = x + + eoi, where x is the true value of 
the CM, Si is a random error distributed as A^(0, sf) and Soi is a systematic error of the Xi 
measurement distributed as A^(0, ctq). Here do is considered as a measure of the scatter of 
the set of systematic errors in input measurements. Evidently, Eqi is unknown, otherwise it 
would be accounted for in the reported value of Xj. We can suppose that Sqi biases Xj but 
does not bias s^. Thus the mathematical expectation of each Xj is 

8{xi) = X + eoi • (9) 

Now we can use the set of n equations ([9]) for i = 1, . . . , n as an equations of condition to be 
solved by the least squares method. As a result of this solution, we obtain an estimate of x 
and its uncertainty ctq, which can be expected close to cr2- Then we can consider cxo ~ o"2 as 
an additive error in the WA uncertainty. Combining this error with ai computed under the 
assumption of absence of systematic errors in Xj we get cXc as defined by Eq [HI 
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Finally, let us notice that (|8]) can be rewritten as 

= J- + ■ (10) 



So, to obtain estimate, there is no need to compute separately both ai and (T2 and then 
use Eq. (E]). 

2.3 Median 



m- 



Another approach routinely used to get the estimate of a CM is computation of a median x. 
The median is known as a robust statistics less influenced by outliers. However its standard 
definition does not provide an estimate of error of a median value (it makes it immune to 
unreliable uncertainties though). 



A possible approach to compute a median uncertainty was proposed by iMiillerl ( 11995 



2000al ). Let Xm be the median of Xi, i.e. Xm = med{xi}. Now we can compute the median 



of the absolute deviations (MAD) as 

MAD = med{\xi- Xm\} . (11) 
The uncertainty of rh is then taken as 

= mad . (12) 

- 1 

One can see that this estimate of the median uncert ainty depends only on the data scatter 



and not on input uncertainties. Later iMiillerl (j2000bl ) proposed a method to take account 



of the uncertainty in input data and thus compute weighted median and its uncertainty. 
However its practical realization, as pointed out by the author, is more cumbersome, and 
the testing results and discussion given therein do not show clear advantage of using weighted 
median. 



3 Tests with artificial data 

In this section, results of two tests with simulated data are presented. These tests were 
constructed to investigate in more details the behavior of the a estimates introduced in 
the previous section. Indeed, many of features discussed here can be seen directly from 
corresponding equations, but not so demonstrative. 

Table[T]shows some numerical examples of computation of WA for two measurements, and 
its standard deviation. To compute we used Q=99%, which corresponds to x^(0.99,l)=6.635 
Using Q=95%, x^(0.95,l)=3.841 does not change main conclusions. Note that, unlike general 
practice, we keep several significant digits in uncertainty just to better show the difference 
between various estimates. 

Classical example 1 shows that (T2 cannot provide reasonable estimate for a whatever 
how large input uncertainty are given. Examples 2-7 and 8-15 show how a estimates change 
with grown Sj for the same Xi. Examples 16-23 show how a estimates change with grown X2 
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for the same Xi and Sj. One can see that no one of di, a2, (J3 provides a satisfactory estimate 
of a for all the examples. 

Several observations from Table [T] are as follows. The ai estimate sometimes it is clearly 
underestimated (examples 8-12). Examples 2-7 and 8-15 illustrate that (J2 (2) cannot pro- 
vide satisfactory estimate, especially in the cases 7, 14, 15, where it seems to be underesti- 
mated. Estimate gives more realistic result, but not in all the cases, e.g. 2-3 and 8-11. 
Moreover, value depends not only on data sample {xj, Si} but also on subjective choice 
of Q. Besides, as can be seen from this test and Eq. [7] that (T3 may show significant jumps 
caused by small changes in input data or confidence level. For these reasons, it was decided 
not to use in further work. 

In contrast to Ui, cr2, and (J3, one can see that cTc approach can provide stable and realistic 
estimate of the standard deviation of the WA. 

In the second test, we use the same set of five measurements with different errors (see 
Fig. [1]). In the upper row the data have minimal uncertainties; in the middle row, all the 
uncertainties are increased by factor of 3; and in the bottom row, all the uncertainties are 
increased again by factor of 3. 

From this Figure, one can see that cxi became greater as grow, as expected. However 
it looks underestimated in case a. The (T2 estimate remains the same for all three examples 
because it does not depend on absolute values of Sj, but only on the dispersion ratio, which 
is the same in all the cases. One would expect however that a should be greater in two 
last examples as compared to the previous ones. Median estimate is close to (T2 as expected 
because both of them depend on the Xi scatter only, and hence am shows the same problems 

as £72. 

The (Jc estimate appears to be optimal because it shows a steady increase from case a to 
case c accounting both for input Si value and Xi scatter. In case a with small Si values o"c is 
determined mainly by the data scatter, and is close to (72. In case b with equal contribution 
of input data uncertainties and scatter (ai = (T2), Uc is just greater by factor of \/2. In case c, 
ac is defined mainly by Sj, which are much greater than the data scatter, and it is close to ai. 
We can say that the ac estimate "automatically" takes account of both input measurements 
scatter and uncertainties without any need in supplement assumptions or parameters like 
significance level. 



4 Application to real data 



In this section, the tests with real data are presented. In the first test, the height differences 
are analyzed between marks 107 and 109 of the local geodetic network of the Svetloe Radio 
Astronomical Observatory of the Institute of Applie d Astronomy, St. Petersburg, Russia 
(iKazarinov and Malkinl . 119971 : iFinkelstein et ai l2006l ) . This analysis includes four levelling 
surveys performed in 1998-2003. Distance between the marks is about 135 m. Results of 
computations are shown in Fig[2j In this case, the scatter of the measurements is rather large 
as compared to measurement uncertainties. For this reason, ai seems to be underestimated. 
Combined estimate as is close to o"2, but may be preferable because accounts also for input 
uncertainties. Median uncertainty is close to a2- In this example, it is difficult to decide 
which of the two latter estimates should be preferred. Both of them looks equally realistic. 
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Table 1: Examples of computation of WA of two measurements xi, X2 with associated 
uncertainties Si, S2- Results of computations are mean x and four estimates of its uncertainty 
c"2, cTs, and cXc computed by (j4]), ([5]), ([7]), and ([8]) respectively. H is computed by ([2]) and 
used to compute (T3 



No. 


Xi 


X2 


S2 


X 


H 


0"! 


0"2 


0"3 


0"c 


1 


1.0 


1.0 


0.5 


1.0 


0.00 


0.354 


0.000 


0.354 


0.354 


2 


1.0 


2.0 


0.1 


1.5 


50.00 


0.071 


0.500 


0.500 


0.505 


3 






0.2 




12.50 


0.141 


0.500 


0.500 


0.520 


4 






0.3 




5.56 


0.212 


0.500 


0.212 


0.543 


5 






0.5 




2.00 


0.354 


0.500 


0.354 


0.612 


6 






1.0 




0.50 


0.707 


0.500 


0.707 


0.866 


7 






2.0 




0.12 


1.414 


0.500 


1.414 


1.500 


8 


10.0 


20.0 


0.1 


15.0 


5000.00 


0.071 


5.000 


5.000 


5.000 


9 






0.5 




200.00 


0.354 


5.000 


5.000 


5.012 


10 






1.0 




50.00 


0.707 


5.000 


5.000 


5.050 


11 






2.0 




12.50 


1.414 


5.000 


5.000 


5.196 


12 






3.0 




5.56 


2.121 


5.000 


2.121 


5.431 


13 






5.0 




2.00 


3.536 


5.000 


3.536 


6.124 


14 






10.0 




0.50 


7.071 


5.000 


7.071 


8.660 


15 






20.0 




0.12 


14.142 


5.000 


14.142 


15.000 


16 


10.0 


10.0 


1.0 


10.0 


0.00 


0.707 


0.000 


0.707 


0.707 


17 


10.0 


11.0 




10.5 


0.50 


0.707 


0.500 


0.707 


0.866 


18 


10.0 


12.0 




11.0 


2.00 


0.707 


1.000 


0.707 


1.225 


19 


10.0 


13.0 




11.5 


4.50 


0.707 


1.500 


0.707 


1.658 


20 


10.0 


14.0 




12.0 


8.00 


0.707 


2.000 


2.000 


2.121 


21 


10.0 


15.0 




12.5 


12.50 


0.707 


2.500 


2.500 


2.598 


22 


10.0 


16.0 




13.0 


18.00 


0.707 


3.000 


3.000 


3.082 


23 


10.0 


17.0 




13.5 


24.50 


0.707 


3.500 


3.500 


3.571 
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a 



xi = 20.0 ± 1.4 
X2 = 18.7 ± 1.7 
X3 = 29.8 ± 1.4 
X4 = 21.1 ± 1.6 
X5 = 23.3 ± 1.0 

Xy,i = 23.00 ±0.60 
x^2 = 23.00 ± 1.81 
x,„r = 23.00 ± 1.91 



Xr 



21.10 ±2.04 



40 
30 

20 t 
10 



1 2 3 4 5 
Measurement number 



xi = 20.0 ±4.2 
X2 = 18.7 ±5.1 
= 29.8 ±4.2 
Xi = 21.1 ±4.8 
X5 = 23.3 ±3.0 

x^i = 23.00 ± 1.81 
x,^2 = 23.00 ± 1.81 
x^c = 23.00 ± 2.56 
= 21.10 ± 2.04 



40 
30 

20 
10 



1 2 3 4 5 
Measurement number 



xi = 20.0 ± 12.6 
X2 = 18.7 ± 15.3 
X3 = 29.8 ± 12.6 
X4 = 21.1 ± 14.4 
X5 = 23.3 ± 9.0 

x^i = 13.00 ±5.42 
x^2 = 13.00 ± 1.81 
x^^ = 13.00 ±5.71 
Xm = 11.10 ±2.04 



40 

30 



12 3 4 5 
Measurement number 



Figure 1: Testing four methods to compute CM with simulated data. The left group of 
points (discs) represents measurements, and the right group of points (circles) represents 
CM estimates. In all three examples, averaged measurements Xi are the same, but their 
uncertainties differ. Uncertainties in example b are three times greater than those in example 
a, and uncertainties in example c are three times greater than those in example b; x^i is 
the WA computed by Eq. ([1]) with the uncertainty ai computed by (jl]). is the same 
WA with the uncertainty cr2 computed by ([5]), x^c is the same WA with the uncertainty 
computed by (IHD, is the median with the uncertainty (Xm computed by ( fT2|) 
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xi = 3847.5 ± 0.4 
X2 = 3847.9 ± 0.4 
X3 = 3846.6 ± 0.5 
Xi = 3848.1 ± 0.2 

x^i = 3847.83 ±0.16 
x^2 = 3847.83 ± 0.26 
x^c = 3847.83 ± 0.31 
Xm = 3847.70 ± 0.32 



t I 



2000 2002 
Year 



Figure 2: Testing four methods to compute a CA with real data: a case of determination of 
the height difference between two geodetic marks. The designations are the same as in Fig[T] 



xi = 15.0 ±0.8 
X2 = 14.4 ± 1.2 
X3 = 11.3 ± 1.1 
X4 = 14.8 ±0.8 
X4 = 14.5 ± 1.5 

x^i = 14.21 ± 0.44 
x^2 = 14.21 ±0.65 
x^c = 14.21 ± 0.79 
x„ = 14.50 ± 0.28 



12 3 4 
Measurement number 



Xi 
X2 
X3 
X4 
X4 



-10.0 ± 1.2 
-12.0 ± 2.8 
-13.9 ±0.9 
-12.4 ±0.6 
-12.0 ±3.0 



-12.42 ±0.45 
-12.42 ±0.59 
-12.42 ±0.74 
-12.00 ±0.37 



12 3 4 
Measurement number 



Figure 3: Testing four methods to compute a CA with real data: a case of determination of 
the Oort constants A (top) and B (bottom). The designations are the same as in Fig[T] 
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In the sec ond test, results of determination of average values of the Oort constants in 



Klackal (120091 ) are revised (see Fig [3]). For comparison, the author's estimates are 14.2 ± 0.5 
for A, and —12.4 ± 0.5 for B, with uncertainties computed as ai and evidently rounded 
up. This value of the WA uncertainty is likely underestimated as compared to input data 
uncertainties and scatter. Median uncertainty also seems underestimated because it accounts 
only for relatively small scatter, and ignores relatively large measurement uncertainties. In 
this example, o"c estimate again looks the most realistic and corresponding to input data. 



5 Conclusion 

Although computation of a CM, in particular WA, is widely used in data analysis, this 
problem has no definite, unambiguity solution yet. In particular, a very important problem 
in most of applications is to obtain a "realistic" estimate for the CM uncertainty. Both 
underestimation and overestimation are equally undesirable. 

In this study, we mostly investigated several basic approaches to compute a WA. Cur- 
rently used methods for computation of the WA uncertainty do not provide satisfactory 
result for many practical tasks. The classical WA uncertainty estimate ai often yields un- 
derestimated value because it is rigorously justified for unbiased Xi only. Another estimate 
(J2, which can be derived from a least squares solution or by appropriate scaling of the input 
uncertainties, does not take account of input uncertainties, but only of the variance ratio. 

In this paper, we propose a new approach to compute the WA uncertainty, combined 
estimate cxc given by (IHl) and ( fTOj) . which is able to account for both uncertainties and scatter 
of input data. We did not propose a rigorous statistical background for this estimate. How- 
ever, it can be shown that in the case when input values are obtained from a normally 
distributed populations, and each value Xi has a normally distributed systematic error, pro- 
posed estimate dc can be derived from a least squares solution. It appears probable that 
this estimate is suitable for practical use until the deviation from normality is very large. 

Several tests with simulated and real measurements have demonstrated that using cXc is 
a simple and effective approach equally suitable for both discrepant and consistent measure- 
ments. It is also important that it provides realistic a estimate even for very small samples 
of 2-3 measurements. 

As to the median approach, it is known as a more robust CM estima t e, but compu tation 



of its uncertainty is not so straightforward. A simple approach by lMiillerl (jl995U2000al ) seems 



to be not satisfactory in some cases, as fo llows from our tests. In this respect, the bootstrap 



method and extended bootstrap method ( iMacMahon et all . |200J) deserve consideration; this 



methods may be too complicated for routine use though. 

In conclusion, one should never forget the following. Despite the method used, the CM 
uncertainty is only a part of real measurement accuracy, namely a Type A uncertainty accord- 
ing to the standard metrological terminology. By definition. Type A uncertainty is computed 
from data using statistical procedures, while a Type B un certainty is ob tained using sup- 



plement information and theoretical considerations (see e.g. iBucherl (|2004J ). Sometimes, the 
type B uncertainty is evaluated from supplement testing data processing, e.g. using subsets 
of the original data set. It can also be a result not of mathematical computation, but con- 
siderations based on the knowledge of the measurement procedures, observational history. 
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and previous experience. Such procedures were u sed e.g. to derive the best estimates of 
some quantities related to Solar System dynamics (jPitjeva EVi . 120091 ). and to compute real- 
istic errors in radi o source positio ns for the second realization of the International Celestial 
Reference Frame (IMa et all . |2009| ) . 
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